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Abstract 

A new class of models, generalizing Asymmetric Exclusion Process for many parallel interacting 
channels, is proposed. We couple the models with boundary reservoirs, study boundary-driven 
phase transitions and show that usually taken hydrodynamic description fails. The adequate 
hydrodynamic limit is then derived. We support our findings with Monte-Carlo simulations of the 
original stochastic system. 
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I. INTRODUCTION 



Many interesting non-equilibrium physical phenomena (shock wave dynamics, boundary 
driven phase transitions, etc.) can be observed already in simplest models such as driven 
diffusing particle systems, recent review of which can be found e.g. in □ D . These systems 
have a unique ability to feel the dynamics of the boundaries thanks to the presence of a flux of 
particles which brings information from boundaries to the bulk. In absence of a flux, indeed, 
boundaries would play marginal role, as in ordinary equilibrium statistical mechanics, but in 
presence of a flux the boundaries may dominate the bulk with the possibility to give rise to 
boundary driven phase transitions. These phenomena were indeed observed in one-species 
driven diffusive models P, U, 0] and in some models containing more than one species of 
particles (many-species models) [6. 7]. In a more general context, the problem was addressed 
in 0| , where it was shown that the stationary state of systems with one species of particles 
with open boundaries in non-equilibrium depends only on the stationary flux. This implies 
that an extremum principle for the stationary flux can be formulated and the properties of 
the stationary state can be obtained directly from the hydrodynamic limit. In this approach 
the stochastic model is mapped into a viscousless conservation law equation which is ill- 
defined (the corresponding Cauchy problem admits multiple solutions). Since the dynamics 
of stochastic system (Markov process) is unique, the problem of regularization, i.e. how to 
single out the physical solution from a set of multiple solutions, arises. As is well known, this 
problem is resolved by adding a phenomenological small dissipative term (of order e) to the 
conservation law equation and then taking the limit e — > 0. For this purpose, a linear viscous 
term of the form e d 2 u/dx 2 is usually considered, although nonlinear terms of the type e f(u) 
d 2 u/dx 2 , with f(u) being an arbitrary convex function, are also possible. For one-species 
systems this regularization procedure is quite robust in the sense that different choices of the 
viscosity term can be shown to lead to the same physical solution. For many-species models, 
however, this is not trivial. In this case the hydrodynamic limit is is described by a system 
of conservation laws (one for each particle species) and the regularization is usually achieved 
by adding a diagonal diffusion matrix of the form e d 2 u l /dx 2 , to the system. In spite of its 
simplicity, there is no evidence that this approach works also for these more complicated 
systems and the regularization problem of multi-species driven diffusive systems is still open. 

The aim of the present Letter is two-fold. From one side, we introduce a new class of 
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multi-species driven diffusive models which generalize the single chain asymmetric exclusion 
process (ASEP). These models possess nice properties such as product measure stationary 
state and particle-hole symmetry. From a physical point view, they describe the motion of 
particles in multi-channel cables with the particles in adjacent cables creating an effective 
barrier potential for the particles in a given cable to flow. From the other side, we study 
the hydrodynamic limit of these models by showing that the conventional regularization of 
the conservation law equations leads to wrong hydrodynamic results. General arguments 
which support this failure will be provided. On the contrary, we show that an alternative 
regularization, obtained directly from the microscopic dynamics, leads to correct results. 
A detailed comparison of the hydrodynamic predictions (both with conventional and alter- 
native regularization) with direct Monte Carlo simulation of the original stochastic model 
confirms these results. 

Although the mapping from hydrodynamic equations to stochastic processes is not unique, 
it is likely that our alternative regularization may work also for other multi-species driven 
diffusive models. 

II. THE MODEL 

To introduce the model, we consider M > 2 discrete chains on which particles can 
hop preferentially in one direction. The hopping between adjacent chains is forbidden but 
particles can move from a site k to an empty site k + 1 on the same chain with rates that 
depend on the particle configuration at adjacent sites of "neighboring" S chains (the results 
will be qualitatively the same for other possibilities of next-neighboring hopping). Notice 
that such a dynamics does not satisfy detailed balance condition meaning that the system is 
far from thermodynamic equilibrium. The hopping rates for different configurations can be 
expressed in terms of a single parameter f3 measuring the strength of the interaction between 
the chains. Let us denote with r n the hopping rate from site k to site k + 1 on the same 
chain, in presence of a total number n of particles in adjacent chain sites. We restrict our 
search to the models having stationary product measures. In this case, calculation shows 
that the rates r n have to satisfy the condition r n — r n _i = const for any n. Moreover, the 
rates r n turn out to depend only on n and not on the particular configurations the n particles 
can assume. With S chains-neighbors, any given pair of consecutive sites k, k + 1 can have 
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Figure 1: Elementary hopping processes happening in the case of M = 3, S = 2, with their rates. 
Four out of 16 possible configurations are shown. The rates (^Q) depend only on the total number 
of particles in the neighboring four sites. 

from n = up to n = 25 particles-neighbors. An example of the hopping process with the 
corresponding rates is shown in Fig. ^for the case 5 = 2. In general, for a chain having 5 
chains-neighbors [sj], the rates of hopping (normalized to tq = 1) can be parametrized as 

r n = l + n((3- l)/25, n = 0, 1, 25. (1) 

Notice that for (3 = 1 the rates become independent on n (i.e. the interchain interaction 
becomes zero) and the system splits into M parallel uncoupled totally asymmetric exclusion 
processes (TASEP) j^, 0]. For (3 < 1, the rates ((H) monotonically decrease with number 
of particles-neighbors n, i.e. the presence of adjacent particles creates an effective barrier 
potential which slows down the particle movement. 

The driven diffusive models described by ((TJ) possess the following properties, i) Product 
measure. One can easily check that the model admits a stationary distribution in the 
class of product measures. This means that the stationary state of the model is spatially 
uncorrelated, both in the longitudinal and in the transverse direction. This property allows 
to obtain a simple analytic expression for the stationary flux (see below). 

ii) Particle - hole symmetry. It is easy to check that an exchange of particles with holes 
plus the substitution (3 1/(3 leaves the model invariant. Since parameter (3 must be is 
in the range < (3 < oo (the rates in are nonnegative), the present property allows to 
further restrict (3 to the range < (3 < 1 only. In the limit of the strongest interaction (3 = 
the hopping of a particle in completely saturated environment becomes impossible r2s = 0. 
Let us consider some specific examples. 

Case 5=1. In this case there is one "neighboring" chain for each given chain. This can 
happen if there are only two chains (M = 2) or if in the model with M consecutive chains, 
the hopping along one chain depends only on the state of the next chain. Due to the property 
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(i), the stationary flux on a given chain for S = 1 can be computed straightforwardly as 

j(p,a) =p(l-p)(l + G9-l)<7) (2) 
where p and a are average densities of particles on a given chain and on a chain " neighboring" 



to it, respectively. The case M = 2 was considered in 10] where the properties of elementary 
excitations were studied. 

Case S = 2. In this case each pair of sites on two neighboring chains have up to 4 particle- 
neighbors (see Fig.EJ. this situation can be realized such a in the geometry of M consecutive 
chains with periodic boundary conditions M + 1 = 1 in the transverse direction, where each 
chain m has two chains- neighbors mil. The hopping rates are then obtained as 

' 1 if n = 0, 

(/3 + 3)/4 if n = l, 

((3 + l)/2 ifn = 2, (3) 

(3/3 + l)/4 ifn = 3, 
k (3 if n = 4. 

The stationary flux on a chain m in a system with average density pk on the chain k can be 
obtained using property (i), as 

jm = Pm(l - An)(l + ^(/? ~ l)(Pm+l + Pm-l)) (4) 



Case S > 2. For arbitrary number of 5 chains-neighbors to a given one, the stationary flux 
is given by 

j m = p m (l -p m ) | li-(/3- 1) Pfc), (5) 

\ neighbors / 

where the sum is taken over the "neighboring" chains to a chain m, and pk is an average 
particle density on chain k. 

From a physical point of view the case S > 2 corresponds to a coaxial cable with many 
fibers, on which particles move, each fiber having S fibers-neighbors surrounding it. A 
real system is always finite, so that boundary conditions, where particles can enter or exit 
the system must be imposed. To model boundaries, we couple model (JTJ) with stationary 
particle reservoirs of densities p™ on the left boundary of chain m, from where particles 
can be injected , and of the densities p^, on the right boundary of it, where they can be 
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extracted. The rates of extraction and injection are obtained from the boundary densities 



and for specific cases were given in 
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To proceed further, we fix the total number of chains, M, the number of chains-neighbors 
of a given chain S, the right and left boundary densities p^, p™, and study bulk stationary 
densities that is, in the limit when time tends to infinity t — ► oo [12], as the interchain 
interaction varies from (3 = 1 (non interacting case) to (3 = (case of maximal interaction). 

It is known that for an Asymmetric Exclusion process, which is a limit of our models 
for (3 = 1, there is a first order non-equilibrium phase transition at = 1 — p™ > |, 
where discontinuous transition from the low density phase p™ tat = p™ for p% = (1 — p™) — 
to a high density phase p™ tat = p'^ for p^f = (1 — p™) + happens. We will choose one 
of the boundary densities p#,p™ in the vicinity of this phase transition point (for (3 = 1), 
different for different m to exclude possible degeneracies, and look at the dependence of 
Pstat(P) f° r different M,S. Alternatively, we shall look at the dynamic properties of the 
stochastic model and the hydrodynamic equation, comparing the time evolution of an initial 
state after certain time. The results of the Monte Carlo calculations can be seen on the 
graphs FigsEISl alongside with the results of numerical integration of corresponding partial 
differential equations arising from stochastic model in hydrodynamic limit. Derivation of 
these equations is given below. 

III. HYDRODYNAMIC LIMIT 

The naive continuum (Eulerian) limit of our stochastic dynamics on the lattice is a system 
of conservation laws 

9P Z M | dj z z = 12 M (6) 

dt dx dx 2 

with a phenomenological vanishing viscosity e — ► (regularization term) on the right hand 
side added. Here p z (x,t) denote coarse-grained particle densities while j z — the particle 
fluxes given by Eq. (jSJ. The regularization term in Eq. (jSJ leads to the correct answer 
for the initial Riemann problem, as compared to the stochastic model (see ^^) but fails to 
describe the reflection from the boundaries (see [l^, Figs|21ISl)- 

An adequate regularization is obtained by averaging exact lattice continuity equations of 
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Figure 2: Comparison of predictions of MonteCarlo (open circles), the hydrodynamic equations 
Eq. ©(line) and Eq. (filled triangles) for stationary densities of particles in the system as a 
function of interchain interaction [5 — 1. Parameters are: M = 2, S = 1, the boundary densities of 
the chains are p\, p\ = 0.4, 0.01 on the left and p^, p\ = 0.2, 0.8 on the right end. 



the stochastic process 



d A 



Jk-l — 3k 



(7) 



for occupation number operators (n^) — > p(x,t) of any given chain, making the lattice 
constant infinitesimally small k, k + 1 — > x + a. For the case (0), the local flux operator jj, 
between the sites k and k + 1 reads (see [10( for a derivation of (JBJ) in the case M = 2, S = 1): 



j fc = n fc (l - fifc+i) 1 + 



/3 — 1 s \ 

p=l J 



(8) 



Here p = 1, ...S are occupation number operators for a site /c of the chains neighboring 
to a given one. We substitute (JHJ) into (|7|). average, factorize and Taylor expand the latter 
with respect to the lattice constant a according to (fik+i) = p(x,t) + + + ^( a3 )' 
etc. Expanding the resulting equation in powers of a and keeping terms up to the second 
order, we obtain the following hydrodynamic equation 

s 

1,2,. ..M, (9) 



d P z {x,t) dr 



d I dp 2 




P =i 



dt dx dx \ dx \ S 

where fP = p p (x,i) are local densities of the chains neighboring to a given chain Z, e 
— — > If— 
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0; £ -> 2ejj- t . Notice that Eq. © is valid for arbitrary M, S and for M = 2, 5 = 1 
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Figure 3: The time evolution of the stochastic model Q with (3 = 0, M = 3. Comparison of 
average density profiles given by Monte-Carlo simulation (points), hydrodynamic equations 
(lines), and naive viscosity approach ©(broken line). The system has evolved for t = 1600 Monte 
Carlo steps from homogeneous initial condition (n?) = p?, matching the left boundary densities 
p\ = 0.2, p\ = 0.01, p\ = 0.1. The corresponding right boundary particle densities are 0.8/0.9/0.1. 
An average over 10 5 histories is taken. One sees that neither the position nor the level of the shock 
are predicted correctly by the naive viscosity approach (the broken line). 

coincides with the hydrodynamic equation derived in Ref. [lo| . 

To check the validity of this approach we have compared direct numerical integrations 
of (J2J) with Monte Carlo simulations of the original stochastic model described by the rates 
in Eq. As a result we obtain that, while the agreement between the Monte Carlo 

predictions and the numerical integration of Eq. Q is excellent, the comparison with the 
standard hydrodynamic limit gives rise to inconsistencies (see discussion below) both in the 
stationary regime and during relaxation towards stationarity. This is shown in Fig|2 for 
system of two and in Fig 01 for three chains with 5 = 1. From these figures it is clear that 
the choice of the viscosity is indeed crucial to obtain correct results, since different choices 
will produce different answers. We checked that Eq. (jHJ) gives results which coincide, up 
to non-universal boundary layer and interface width, with Monte Carlo simulations also for 
other values of M, S, and we believe this fact is valid in general. 

The above analysis can be easily generalized to the case in which the hopping along the 
chains occurs in both directions. In this case the corresponding flux to the left is obtained 
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from the one to the right by exchanging p(l — p) in (J3J) with the term (p — q)p(l — p), where 
p and q denote the particle hopping rates on the m-th chain to the right and to the left in 
he empty environment (the general conditions for a product measure state can be found in 

3)- 

It is interesting to discuss why a diagonal diffusion matrix leads to inconsistent results for 
multi-species systems. To this regard, we remark that one can formally obtain a diagonal 
diffusion matrix in Taylor expansion of (J7J) by adding to the flux (jHJ) an additional term as, 
e.g. 

* = Mi - <w (i + ^± K + *w) + (>w, - <w £ K + *U 

(10) 

If one substitutes the above expression into (J7J) and performs the same analysis as before, 
one obtains the resulting equation in the form (jHJ). In this case, however, the factorization 
becomes invalid because the flux operator (jl(J|) describes different process with correlations in 
the stationary state and therefore with a stationary flux different from Thus, by making 
the diffusion matrix diagonal one expects inconsistencies. This can be seen by considering a 
microscopic state with all but one chains completely filled with particles, the chain partially 
filled having an empty region separated by a completely filled one, = 0, k < and = 1, 
k > 0. According to the dynamic rules in Figd no movement is allowed, while from Eq. 
(J7J) (with averaged and factorized flux as in Eq. (jTU|l ). we obtain that ^n ^ 0. Notice that 
although the choice of the second term in Eq. (fTU|) is not unique, the above arguments for 
inconsistency would still be valid in presence of alternative choices. These inconsistencies 
are clearly seen on Fig0 where a comparison between stationary densities of the two chain 
system (M = 2) obtained from Monte Carlo calculations and hydrodynamic equations ()6I9|) 
is made, as a function of the interchain interaction 1 — (3. The naive description © shows a 
phase transition at (5 ~ 0.1 , while Monte Carlo and © show no effect. The occurrence of 
a phase transition to high densities can be understood as follows. Since the right boundary 
density of particles of chain 2 is high {{p2) R = 0.8, it becomes more and more difficult for 
the particles of chain 1 to exit the chain with increasing interaction, so that "traffic jam" 
occurs at the exit, this leading to the phase transition. This effect is greatly overestimated 
by the naive hydrodynamic description ( EJ) (in physical system the transition happens for 
higher left boundary density (pi) L > 0.48). 
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IV. CONCLUSION 



In conclusion, we have introduced a class of models which generalize the asymmetric 
simple exclusion process for an arbitrary number of chains and with hopping rates given as 
a function of the local configuration in the neighboring chains (J!Q) describing the effective 
friction a particle encounters moving in a dense environment. For these models we derived 
the stationary flux, and modified the conventional choice of the viscosity term in the hydro- 
dynamic limit (see, e.g., 15]) in order to get a good agreement with Monte Carlo simulations 
of the original stochastic process. Although the problem of the hydrodynamic limit has been 
discussed on the specific example of the introduced models, the results are expected to be 
valid in general (the advantage of using our model is only for computational convenience 
since the product measure property allows to obtain the bulk flux, boundary conditions, 
etc., in an easy and straightforward manner.) 
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